1つ上の潜在成長曲線モデル

竹林 由武(統数研)

2016-3-27 @Hirosima.Univ

Topics

 

1. 非線形の成長曲線モデル

  • polynominal model
  • latetn basis model
  • piecewise models

 

2.多変量の成長曲線モデル

  • parallel process growth model
  • autoregressive latent trajectory Model
  • latent change score model

非線形モデルのモチベ

2次(以上)の曲線で

polynominal model

データドリブンに

latent basis model

ある時点から傾きが異なる

piecewise model

サンプルデータ

ポジ感情の経時データ

X1 X2 X3 X4 X5 X6
mean 3.49 3.30 3.27 3.38 3.35 3.30
sd 0.97 0.98 0.97 0.98 0.98 0.99
time 0.00 1.00 2.00 3.00 4.00 5.00

Garland, E. L., Geschwind, N., Peeters, F., & Wichers, M. (2015). Mindfulness training promotes upward spirals of positive affect and cognition: multilevel and autoregressive latent trajectory modeling analyses. Frontiers in psychology, 6.

Rで潜在成長曲線モデル

lavaanでまじめに書く

level =~ 1* bmi1 +1* bmi2 +1* bmi3 +
         1* bmi4 +1* bmi5 +1* bmi6 
slope =~ 0 * bmi1 + 1 * bmi2 + 2 * bmi3 + 
         3 * bmi4 + 4 * bmi5 + 5 * bmi6 
bmi1 ~~(vare)* bmi1 
bmi2 ~~(vare)* bmi2 
bmi3 ~~(vare)* bmi3 
bmi4 ~~(vare)* bmi4 
bmi5 ~~(vare)* bmi5 
bmi6 ~~(vare)* bmi6 

結構めんどい…

救世主登場

RAMpathパッケージ

  • 潜在成長曲線系のlavaanコードを自動生成して実行してくれる関数が充実
    • 成長曲線モデル: ramLCM
    • 潜在差得点モデル:ramLCS
    • 2変量の潜在差得点モデル: ramBLCS
  • lavaanのモデルを吐き出せるので、RAMpathで実行した後lavaanで微調整という流れで効率よくモデリング

ramLCM

  • 切片のみのモデル (model=‘no’)

  • 線形モデル (model=‘linear’)

  • 二次曲線モデル (model=‘quadratic’)

  • latent basisモデル (model = ‘latent’)

一行で簡潔に推定

一気にどん

library(RAMpath)
fit.all<-ramLCM(data=data,outcome=1:6, model='all')

個別にどん

fit.no<-ramLCM(data=data,outcome=1:6, model='no')
fit.linear<-ramLCM(data=data,outcome=1:6, model='linear')
fit.quadratic<-ramLCM(data=data,outcome=1:6, model='quadratic')
fit.latent<-ramLCM(data=data,outcome=1:6, model='latent')

切片のみモデル

lavaanコード

cat(fit.all$model$no)
## level =~ 1* X1 +1* X2 +1* X3 +1* X4 +1* X5 +1* X6 
##  X1 ~~(vare)* X1 
##  X2 ~~(vare)* X2 
##  X3 ~~(vare)* X3 
##  X4 ~~(vare)* X4 
##  X5 ~~(vare)* X5 
##  X6 ~~(vare)* X6
  • 傾き因子を排除
  • 初期値から得点の変化なし

線形モデル

lavaanコード

cat(fit.all$model$linear)
## level =~ 1* X1 +1* X2 +1* X3 +1* X4 +1* X5 +1* X6 
##  slope =~  0 * X1 + 1 * X2 + 2 * X3 + 3 * X4 + 4 * X5 + 5 * X6 
##  X1 ~~(vare)* X1 
##  X2 ~~(vare)* X2 
##  X3 ~~(vare)* X3 
##  X4 ~~(vare)* X4 
##  X5 ~~(vare)* X5 
##  X6 ~~(vare)* X6
  • 傾き因子の因子負荷が単調増加 (0, 1, 2, 3, 4, 5)

二次曲線モデル

lavaanコード

cat(fit.all$model$quadratic)
## level =~ 1* X1 +1* X2 +1* X3 +1* X4 +1* X5 +1* X6 
##  slope =~  0 * X1 + 1 * X2 + 2 * X3 + 3 * X4 + 4 * X5 + 5 * X6 
##  quadratic =~  0 * X1 + 1 * X2 + 4 * X3 + 9 * X4 + 16 * X5 + 25 * X6 
##  X1 ~~(vare)* X1 
##  X2 ~~(vare)* X2 
##  X3 ~~(vare)* X3 
##  X4 ~~(vare)* X4 
##  X5 ~~(vare)* X5 
##  X6 ~~(vare)* X6
  • 二つ目の傾き因子の負荷は, 一つ目の傾き因子の二乗 (0, 1, 4, 9, 16, 25)

latent basisモデル

lavaanコード

cat(fit.all$model$latent)
## level =~ 1* X1 +1* X2 +1* X3 +1* X4 +1* X5 +1* X6 
##  slope =~  0 * X1 +start( 1 )* X2 +start( 2 )* X3 +start( 3 )* X4 +start( 4 )* X5 + 5 * X6 
##  X1 ~~(vare)* X1 
##  X2 ~~(vare)* X2 
##  X3 ~~(vare)* X3 
##  X4 ~~(vare)* X4 
##  X5 ~~(vare)* X5 
##  X6 ~~(vare)* X6
  • 傾きの因子負荷を自由推定する
  • データドリブンな遷移パターン -ramLCMの自動コードは若干ミスってる

コードミスを修正

modelXlat<-'
level =~ 1* X1 +1* X2 +1* X3 +1* X4 +1* X5 +1* X6 
slope =~  0 * X1 +start( 1 )* X2 +start( 2 )* X3 +start( 3 )* X4 +start( 4 )* X5 + 5 * X6 
X1 ~~(vare)* X1 
X2 ~~(vare)* X2 
X3 ~~(vare)* X3 
X4 ~~(vare)* X4 
X5 ~~(vare)* X5 
X6 ~~(vare)* X6 
'
Xlat<-lavaan::growth(modelXlat,data)
Xlat.fit<-round(fitmeasures(Xlat)[c("chisq","df","pvalue","cfi","srmr","rmsea","aic","bic")],digits=2)

遷移プロット

  • plot.growth関数(自作)
source("script/plot.growth.R")
a<-plot.growth(fit.all, type="no")+theme_bw()
b<-plot.growth(fit.all, type="lin")+theme_bw()
c<-plot.growth(fit.all, type="quad")+theme_bw()
d<-plot.growth(fit.all, type="latent")+theme_bw()

適合度

fits<-round(fit.all$fit[
            c("chisq","df","pvalue","cfi",
              "srmr","rmsea","aic","bic"),],digits=2)
fits[,2]<-Xlat.fit
datatable(fits,option=list(dom='t'))

ポチるおまけ

推定結果の出力

parm<-parameterEstimates(fit.all$lavaan$quadratic)
parm[,5:10]<-round(parm[,5:10],digits=3)
datatable(parm[c(28:30,37:39),],
          extensions = 'Scroller', 
          options = list(dom= ' t',
                         deferRender = TRUE, 
                         scrollY = 200, 
                         scroller = TRUE))

推定平均値の遷移プロット

source("script/plot.growth.R")
g<-plot.growth(fit.all, type="quad")
g+ylim(3,4)+theme_bw()

Piecewise Model

 コード全体

model1 <-'

#切片因子の設定
i =~ 1*X1 + 1*X2 + 1*X3 + 
1*X4 + 1*X5 + 1*X6

#傾き因子の設定
s1 =~ 0*X1 + 1*X2 + 2*X3 +
3*X4 + 3*X5 + 3*X6

s2 =~ 0*X1 + 0*X2 + 0*X3 +
0*X4 + 1*X5 + 2*X6


#切片と傾きの分散
i ~~ i ; s1 ~~ s1 ; s2 ~~ s2; 

#因子間相関
i ~~ s1 + s2; s1 ~~ s2

#因子平均
i ~ 1 ; s1 ~ 1 ; s2 ~ 1

#誤差分散
X1 ~ 0; X2 ~ 0; X3 ~ 0
X4 ~ 0; X5 ~ 0; X6 ~ 0
'

切片因子の設定

  • 因子負荷を1に固定
#lavaan code
i =~ 1*t1+1*t2+1*t3+1*t4+1*t5

前半の傾きの設定

前半の傾き(s1)の因子負荷を
区分時点以降同値に固定

#lavaan model code
i=~0*t1+1*t2+2*t3+3*t4+3*t5+3*t6

後半の傾きの設定

後半の傾き(s1)の因子負荷を
区分時点まで0に固定

#lavaan model code
i=~0*t1+0*t2+0*t3+0*t4+1*t5+2*t6

その他の設定

  • 切片と傾きの分散を自由推定
  • 因子間相関を自由推定
  • 因子平均を自由推定
  • 誤差分散を0に固定
#切片と傾きの分散
i ~~ i ; s1 ~~ s1 ; s2 ~~ s2
#因子間相関
i ~~ s1 + s2 ; s1 ~~ s2
#因子平均
i ~ 1 ; s1 ~ 1 ; s2 ~ 1
#誤差分散
bmi1 ~ 0; bmi2 ~ 0; bmi3 ~ 0
bmi4 ~ 0; bmi5 ~ 0; bmi6 ~ 0

推定の実行

library(lavaan)
model1.fit<-lavaan::growth(model1, data)

フィットしてる?

fit1.m<-round(fitMeasures(model1)[c("chisq","df","pvalue",
"cfi","srmr","rmsea")],digits=2)
fit1.m<-t(as.data.frame(fit1))
print(xtable(fit1.m),comment=F,type="html")
  • piecewise

  • quadratic

前半と後半の傾きに違いがあるか?

  • 傾き因子の平均に等値制約を置いたモデルと比較

    #傾き因子平均が等値
    s1 ~ (a)*1 ; s2 ~ (a)*1
# model1.2fit = 等値制約のモデル
anova(model1.2fit, model1.fit)
## Chi Square Difference Test
## 
##             Df   AIC   BIC  Chisq Chisq diff Df diff Pr(>Chisq)
## model1.fit  12 26252 26336 324.30                              
## model1.2fit 13 26251 26330 325.46     1.1551       1     0.2825

もうちょい頑張る

  • 前半を2次曲線モデルに
model2 <-'
#切片因子の設定
i =~ 1*X1 + 1*X2 + 1*X3 + 1*X4 + 1*X5 + 1*X6

#傾き因子の設定
s1 =~ 0*X1 + 1*X2 + 2*X3 +3*X4 + 3*X5 + 3*X6
s2 =~ 0*X1 + 1*X2 + 4*X3 +9*X4 + 9*X5 + 9*X6
s3 =~ 0*X1 + 0*X2 + 0*X3 +0*X4 + 1*X5 + 2*X6

#切片と傾きの分散
i ~~ i ; s1 ~~ s1 ; s2 ~~ s2; s3 ~~ s3; 

#因子間相関
i ~~ s1 + s2 + s3; s1 ~~ s2 + s3 ; s2 ~~ s3 ;

#因子平均
i ~ 1 ; s1 ~ 1 ; s2 ~ 1 ; s3 ~ 1

#誤差分散
X1 ~ 0; X2 ~ 0; X3 ~ 0; X4 ~ 0; X5 ~ 0; X6 ~ 0
'

 どやさ

-今回の適合度

model2.fit<-lavaan::growth(model2, data=data)
fit2.m<-round(fitMeasures(model2.fit)
[c("chisq","df","pvalue",
"cfi","srmr","rmsea")],digits=2)
fit2.m<-t(as.data.frame(fit2.m))
datatable(fit2.m,options=list(dom="t"))

  • さっきのpiecewiseモデル

  • 圧倒的なfit感!!

圧倒的なfit感!!

 期間ごとに切片が異なるモデル

  • 二つ切片因子を設定
model2.2 <-'
i1 =~ 1*X1 + 1*X2 + 1*X3 + 0*X4 + 0*X5 + 0*X6
i2 =~ 0*X1 + 0*X2 + 0*X3 + 1*X4 + 1*X5 + 1*X6
s1 =~ 0*X1 + 1*X2 + 2*X3 + 3*X4 + 3*X5 + 3*X6
s2 =~ 0*X1 + 0*X2 + 0*X3 + 0*X4 + 1*X5 + 2*X6

i1 ~~ i1
i2 ~~ i2
s1 ~~ s1
s2 ~~ s2

i1 ~~ i2 + s1 + s2
i2 ~~ s1 + s2 
s1 ~~ s2
i1 ~ 1
i2 ~ 1
s1 ~ 1
s2 ~ 1

X1 ~ 0
X2 ~ 0
X3 ~ 0
X4 ~ 0
X5 ~ 0
X6 ~ 0 '

こんな感じでプロット

成長曲線モデルまとめ

  • 遷移パターンは複数のモデルを比較検討
  • 切片のみ、線形、二次曲線、latent basis
  • RAMpath パッケージで瞬殺
  • 区切りが明確な場合はpiecewiseも試してみよう

2変量の関係性が知りたい

2つの時系列データ

  • ポジ感情とポジ認知の経時データ

Garland, E. L., Geschwind, N., Peeters, F., & Wichers, M. (2015). Mindfulness training promotes upward spirals of positive affect and cognition: multilevel and autoregressive latent trajectory modeling analyses. Frontiers in psychology, 6.

有効なモデル

  • parallel process growth model
  • auto regressive latent trajectory model
  • latent change score model

parallel process growth model

 

  1. 各変数に最適な潜在成長曲線モデルを検討
  2. 潜在変数(切片と傾き)間の相関を推定

ALTM

ALTM

2つの時系列データ

交差遅延効果モデル

成長曲線モデル

交差遅延効果+潜在成長曲線=ALTM

自己回帰のパス

-変数の安定性

交差遅延のパス

  • 変数間の因果関係

 切片と傾き (全体の遷移パターン)

LCS